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FULL-SCALE DIRECT NUMERICAL SIMULATION OF 
TWO- AND THREE-DIMENSIONAL INSTABILITIES AND RIVULET FORMATION 

IN HEATED FALLING FILMS 1 
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+ Department of Mechanical Engineering and Materials Science 
Rice University, Houston, Texas 77251-1892, U.S.A. 

* Department of Mechanical Engineering 
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A thin film draining on an inclined plate has been studied numerically using finite element 
method. Three-dimensional governing equations of continuity, momentum and energy with 
a moving boundary axe integrated in an Arbitrary Lagrangian Eulerian frame of reference. 
Kinematic equation is solved to precisely update interface location. Rivulet formation based 
on instability mechanism has been simulated using full-scale computation for the first time 
in the literature. Comparisons with long-wave theory are made to validate the numerical 
scheme. Detailed analysis of two- and three-dimensional nonlinear wave formation and 
spontaneous rupture forming rivulets under the influence of combined thermocapillary and 
surface- wave instabilities is performed. 
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Long- wave Theory 
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1. INTRODUCTION 


The flow of liquid-film on a solid substrate has many significant engineering applica- 
tions in material processing, biomedical engineering, and nuclear, aerospace and chemical 
industries. Some of the technological importance are tertiary oil recovery; processes as- 
sociated with multi-phase flow through porous media; spreading of liquid as in coating 
process; coalescence of drops and bubbles in foams and emulsion; fabrication of chips in 
micro-electronics; study of cancer cells; development of anti-icing system for aircraft wings. 

The most widely observed phenomena in thin-film flows are those caused by the inter- 
facial instability such as formation of transverse roll waves on the surface, longitudinal roll 
patterns, wave breaking, rupture, breaking of stream of liquid into independent rivulets, 
evaporation and termination of liquid layer at a contact line and formation of dry spots 
Various hydrodynamics forces viz. hydrostatic pressure, surface-tension, inertia, thermo- 
capillary and vapor recoil compete with each other and consequently decide the stability of 
the system. These instabilities can also be described quantitatively by the film thickness, 
amount of heating, heat loss at the interface and the angle of inclination of the plane. We 
follow Goussis and Kelly (1991) and define three modes of instabilities viz P-mode, S-mode 
and H-mode that commonly occur m this type of flow and neglect evaporation, surface con- 
tamination and inter-molecular forces. When the plate is tilted, the liquid drains down due 
to gravity. If the inertia of the mean flow dominates the stabilizing effects of hydrostatic 
pressure, hydrodynamic or H mode instability sets m. This was identified by Yih (1955) 
and Benjamin (1957) who formulated linear stability problem by analyzing the growth of 
spatially periodic two-dimensional disturbances in laminar free-surface flow They arrive 
at the following condition for neutral stability: 


gd 0 3 sinl3 

2u 2 



(1) 


where, g is the gravitational acceleration, do is the mean film thickness, u is the kinematic 


2 



viscosity of the fluid and fi is the angle of inclination. In the limit, as (3 — * 7r/2, d 0 — >0, i e. 
the flow is always unstable. 

When the plate is horizontal and is heated, the interface can be subjected to temper- 
ature gradient in any arbitrary orientation. There exists a purely a static base state if the 
gradient is imposed in the normal direction. Instabilities of this state was first studied by 
Pearson (1958) and is called as P-mode. If the gradients are imposed along the interface, 
no static state exists and the thermocapillary can drive a steady shear flow whose insta- 
bilities are often times periodic hydrothermal waves (Smith and Davis, 1983a, b). P-mode 
occurs when the convection is significant and the wavelength of the imposed disturbance is 
comparable to mean film thickness. The following condition must be satisfied for instability: 


T^= > 32.073 
al ph 


(2) 


where a is the surface-tension of the liquid, T is temperature, Cp is the heat capacity, p 
is the density, p is the dynamic viscosity and h is the heat-transfer coefficient. For this 
mode of instability to develop, the energy transferred to the disturbance and the work done 
by thermocapillary forces have to be large enough to overcome the losses due to viscous 
dissipation and surface heat transfer. This mode does not require surface deformation. 

S-mode is solely induced by surface deformation when the thickness of the film is small. 
This was first studied by Scnven and Sterlin (1964). Here, surface-tension suppresses 
disturbances of shorter wavelength and the film becomes unstable only to long surface 
waves. The balance between thermocapillary force and stabilizing hydrostatic forces can 
be expressed as: 



( 3 ) 


where k is the thermal conductivity of the liquid. As summarized by Goussis and Kelly 
(1991), at sufficiently thin layers the hydrostatic pressure dominates the effect of inertia 
and the film is stable with respect to H-mode. As the depth of the layer increases the 
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inertia dominates and the film becomes unstable. Similarly, for thm layer, the interaction 
of the basic temperature with perturbation velocity can not transfer enough energy to cause, 
instability. However, in thick fluid layer such transfer of energy is possible and the film is 
unstable to P-mode. On the other hand, when the layer is thin, thermocapillary dominates 
hydrostatic pressure and the film is unstable with respect to S-mode However, as the layer 
thickness increases, due to hydrostatic stabilization, film may become stable In this study, 
we focus only on S-mode type instability and neglect surface contamination, evaporation 
and inter-molecular forces. 

The presence of both thermocapillary and surface- wave instabilities can cause a heated 
falling film to rupture and eventually form rivulets. A rivulet is a stream of liquid flow- 
ing down a solid surface and sharing an interface with a surrounding gas (Young and 
Davis, 1987). Hence, it is very important to understand the mechanism of spontaneous 
rupture before we study the dynamics of rivulet formation. In this regard, we first ana- 
lyze spontaneous rupture in truly two-dimensional flow and then extend the analysis to 
three-dimensional rivulet formation. 

Theoretically, the combined thermocapillary and surface-wave instability can be stud- 
ied using linear or weakly nonlinear analysis. Linear stability analysis on these type of flows 
has been performed since the work of Lm (1975). Kelly et al. (1986) identified a “stability 
window” below and above which the flow becomes unstable due to thermocapillary and 
surface- wave instability respectively. The window exists due to the stabilizing effects of 
hydrostatic pressure. Goussis and Kelly (1991) also showed that these instabilities can 
reinforce each other and a disturbance takes the form of a transverse wave when the film 
is very thin and a longitudinal roll wave when it is moderately thick. Linear analysis also 
gives useful information on critical layer thickness, inclination angle, amount of heating 
and cut off wavenumber for neutral stability. However, we can not follow the dynamics of 
the flow which is the focus of the present study. 
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Since the instabilities appear in the form of long interfacial waves, long-wave evolu- 
tion equation of Benny (1966) type are very useful Burelbach et al. (1988) considered 
sufficiently thin horizontal layer and studied long- wave instabilities in the presence of evap- 
oration, vapor recoil and van der Waals forces. Joo et al. (1991) generalized this study 
to include the effect of mean flow in the absence of van der Waal forces. They studied 
the nonlinear flow development by numerically integrating the evolution equation. They 
followed the flow up to the point of rupture when thermocapillary is significant and up to 
the point of wave-breaking when surface- wave instability is dominant and showed that the 
rupture always follows a characteristic “fingering” process and substantial local thinning. 
This process is very sensitive to initial condition. Later on, Joo and Davis (1992) extended 
the analysis to three-dimensional isothermal flows on a vertical plane and identified a new 
secondary instability in which three-dimensional disturbance is spatially synchronous with 
two-dimensional wave. The instability grows for sufficiently small cross-stream wavenum- 
bers and does not require any threshold amplitude. In addition, they studied the three- 
dimension layers by posing various initial value problems and numerically integrating the 
long-wave evolution equation. Recently, Joo et al (1995) included thermocapillary also in 
their analysis and demonstrated a mechanism of rivulet formation solely based on instabil- 
ity phenomenon. They showed that the film first ruptures by “fingering” mechanism and 
then forms rivulets and when thermocapillary and surface-wave instabilities are properly 
balanced different flow patterns can be observed. Though long-wave evolution equation 
can predict the permanent wave form behavior and follow the evolution of finite amplitude 
disturbances, toward rupture the inertial forces assumed small become significant and the 
slope of the interface increases continuously. These phenomena eventually violate the basic 
lubrication type approximation used m long-wave theory formulation. Hence, to study the 
complicated nonlinear flow development without any a prion assumptions, the complete 
system of Navier-Stokes equations must be solved. 
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Due to the irregular and time varying domain involved, Finite Element Method is 
the most popular choice as a numerical tool to investigate thm film flows So far most 
of the numerical simulations done focus on isothermal flow (Bach and Villadsen 1984; 
Kheshgi and Scriven, 1987, Ho and Patera 1990; Malamataris and Papanastasiou 1991; 
Salamon et al. 1994; Chippada 1995). Recently Krishnamoorthy et al. (1995) have studied 
rupture dynamics in two-dimensional non-isothermal flows. Bach and Villadsen (1984), 
Kheshgi and Scriven (1987) and Malamataris and Papanastasiou (1991) use a Lagrangian 
Finite Element Method to handle moving boundary and control mesh distortion through 
rezoning. Kheshgi and Scriven (1987) used Galerkin weighted residual, implicit predictor 
corrector, mixed finite element formulation and studied isothermal thm film flows. Ho and 
Patera (1990) studied the stability of these flows using Legendre spectral element method. 
They used Orr-Sommerfeld theory and experimental studies of Kapitza and Kapitza (1949) 
to compare their spectral element calculations. Salamon et al. (1994) used finite element 
equations written m a reference frame translating at wave speed to study finite amplitude 
waves propagating at constant speed. They found good agreement with long-wave theory 
for small amplitude waves, but found their results to qualitatively diverge from long-wave 
results for large amplitude waves. They also studied the nonlinear interaction between the 
waves and the secondary subharmonics bifurcation to longer waves. Krishnamoorthy et al. 
(1995) solved the governing equations written in Arbitrary Lagrangian Eulerian frame of 
reference and showed that m two-dimensional heated film, “fingering” process leading to 
rupture is not an artifact of long-wave theory but it is an actual phenomena. 

In the present study, we the solve complete system of governing equations along 
with suitable boundary conditions m an ALE frame of reference for both two- and three- 
dimensional thm-film flows and study the nonlinear flow development to spontaneous rup- 
ture and rivulet formation. First, mathematical formulation is discussed in Section (2). 
The numerical scheme is explained in Section(3). In section(4) results of our simulation 
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are discussed and comparison with long-wave theory is made wherever possible. Some 
of the issues that arise during the course of this analysis are discussed and concluded in 
Section(5). 

2. MATHEMATICAL FORMULATION 
2.1 Arbitrary Lagrangian Eulerian (ALE) Formulation 

Thm film flow problem involves more comprehensive analysis of the local flow pattern 
and we need to solve free-surface Navier-Stokes equations in primitive variable formulation. 
We can not impose any ad hoc restriction on the pressure distribution. This results in 
additional nonlinearities associated with the geometry of the free-surface. To solve the 
problem more efficiently, we use an arbitrary Lagrangian Eulerian (ALE) frame of reference 
and write the governing equations accordingly. 

The ALE description of the fluid flow is called referential kinematic description of the 
flow, since the governing equations are written in a frame of reference that move indepen- 
dent of the fluid motion. This formulation was initially proposed by Hirt et al (1974) 
and later on used by many researchers (Chan, 1975, Hughes, Liu and Zimmerman 1981, 
Ramaswamy and Kawahara 1987, Ramaswamy 1990, Soulaimani et al. 1991, Chippada et 
al. 1995) in modeling free-surface flow problems and fluid-structure interaction problems 
(Donea et al. 1982, Donea 1983, Liu et al. 1988). ALE formulation has been derived by 
Donea et al. 1982, Ramaswamy and Kawahara (1987), Soulaimani (1991), Lacroix and 
Garon (1992), among many others and is briefly described next (Chippada 1995). 

Let B 0 be the open region occupied by fluid particles at t=0 as shown in Fig.l. This 
is also called material domain. The position vector of a point P in Bo is denoted by 
X t =(Xi : X 2 ,Xo). Bt is the open region occupied by Bq after some time f>0. The point P 
occupies a unique point p m B t whose position vector is denoted by x t ={x\^X 2 ^xz). It is 
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material domain 


spatial domain 



referential domain 


Figure 1: Arbitrary Lagrangian Eulerian description of the flow. Figure adopted from 
Chippada (1995) 

assumed that the mapping between P and p is continuous, unique and invertible 


x t = X t = <f> 


(4) 


Lagrangian (material) velocity is defined as 


d 


«*„()= 


(5) 


and Lagrangian acceleration is defined as 

o2 

4>{x t ,t) = — 4>(x t ,t) 

In terms of spatial coordinate x t , the Eulerian velocity is defined as 

Uj(x t ,t) = ^ <f>{x tl t) 

and the Eulerian acceleration can be shown to take the form: 

®t(®xj^) = “b U]{x , , t )ti, j (Xj , t) 


( 6 ) 

(7) 

(8) 
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Unlike pure Eulerian description, pure Lagrangian description does not involve any nonlin- 
ear convective terms. However, Eulerian description is most widely used in fluid mechanics 
problems 

In the ALE formulation, a third domain called the referential domain is specified. 
At time t— 0, it occupies an open region i? 0 - This body has its own motion and at some 
time later, t> 0, it coincides with the material body which has moved to B t . That is, the 
referential point q whose position vector is x,=(xi,X 2 , x 3 ) coincides with point p in space 
coordinate after some time t> 0. The mapping between these two regions is given by 

x t = \ (x t ,t) (9) 

and the referential velocity ( g t ) will be 

= ^X(x t1 t) ( 10 ) 

The relative motion of the fluid with respect to the referential frame is expressed as 

x, = A _1 (x,,t) = A - l {<f>(X t ,t),t) = (11) 


and the relative velocity is 




X , 


(12) 


The material point P and the referential point q arrive at the spatial point p through 
independent motions which are related as follows - 


x t = X(x t ,t) = A (ip(X t ,t),t) = (f>(X t ,t) 


(13) 


From the above relation, we can derive 


(®»J^) ^t(*^»)t) 


Xx 

g,(x t: t) + F tJ (x t , t)V?(x t , t ) 


(14) 
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(15) 

(16) 



where, F t;) is called gradient deformation tensor defined as 


= (17) 

The spatial acceleration m the referential coordinates can be shown to be 

= u,,t + V/u t ,j(x„i) (18) 

= u*,t + F-\iij - gj)u t , 3 (x t ,t) (19) 

Comparing spatial acceleration (Eq.8) with referential acceleration (Eq.19), the difference is 

that gradients are with respect to referential coordinate and the spatial velocity is replaced 
by the relative velocity. Equation (19) can also be interpreted as material conservation 
laws with respect to arbitrary moving points. In the event a grid point coincides with 
the material point, the relative velocity becomes zero Consequently, the set of 

equations become Lagrangian. Similarly, a pure Eulerian description can be obtained by 
setting g t to zero. The ALE approach combines the advantages of both Lagrangian and 
Eulerian methods and avoids mesh distortion. In our problem, the referential motion is 
related to the fluid motion and at the free boundary, the mesh points are moved normal 
to the interface with fluid velocity to prevent the loss or gain of fluid material. We use a 
time stepping procedure in which referential velocity computed from previous time step is 
used. With this simplification, we write the governing equations. 

2.2 Governing Equations 

Consider a thin film of Newtonian, incompressible, non-volatile, constant property 
(density p, viscosity p, thermal conductivity k , thermal diffusivity a) liquid kept on a plate 
maintained at a constant temperature T w and inclined at an angle /3 to the streamwise 
direction. The film is thick enough so that the continuum theory is valid and neglect the 
buoyancy forces. It is unbounded in both streamwise and spanwise directions, but bounded 
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Vapor 



Figure 2: Physical configuration of a thin film flowing down a heated inclined plate (two- 
dimensional sectional view). 

above by a passive gas (having negligible density and viscosity) of far-field temperature 
T s (< T w ) and zero pressure. A two-dimensional cross-sectional view of the physical domain 
is shown in Fig. 2. The heat conducted across the liquid layer is lost through the interface 
due to the convection and affects the interfacial temperature T t . We assume that surface 
tension decreases monotonically with temperature: 

<T) = [1 - 7 (T - T s )] , (20) 

where <To is the value of the surface tension at the reference temperature T s and 7 is 
the temperature coefficient. The thermo capillary induced by surface tension gradient is 
measured by Marangoni number M—^(T W — T s )do/2ga. In an undisturbed layer of liquid 
film, the streamwise component of the velocity reaches maximum at the interface and is 
expressed as gd 0 2 s'm/3/u where g is gravitational constant, d 0 is the initial mean thickness 
of the film and u = p/ p is the kinematic viscosity of the liquid. The Reynolds number 
based on the initial maximum velocity and the mean thickness is G sin/?, where G=gdo 3 / v. 
This is also called Galileo number in some literature. This parameter is a measure of the 
film thickness. 
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The governing equations are time- dependent and three-dimensional conservation laws 
for mass, momentum and energy Using mean film thickness ( d Q ) and viscous time (do 2 /v) 
as the scales of motion, the non-dimensional form of the governing equations in an ALE 
frame of reference is written as: 


U M = 0, 

(21) 

Q t + ( u j g j) u i,3 - + G5 t 3 

(22) 

d0 1 

-qI + K - &) e ,3 = 

(23) 


Here, u,=(u,u,i£j) and g t ={g x ,g y ,gz) are respectively the velocity vector and the grid-pomt 
velocity vector, P(=z//a) is the Prandtl number, 9=(T —T s ) / (T w —T s ) is the non-dimensional 
temperature, i= 1,2,3 andj=l,2,3 , denotes the partial derivative and 8 tJ is the Kronecker 
delta. The stress tensor cr tJ is expressed as: 

<r tJ = -p8 tJ + (u t)J + u J>t ) (24) 

where p is pressure. The above equations refer to a right-handed Cartesian coordinate 
system x,=(x,y,z) whose origin is on the plate, z-axis is aligned with streamwise direction, 
^-axis runs along spanwise direction and z-axis points normal to the plate into the liquid. 

The inclined plate is maintained at a constant temperature T w and no-slip boundary 
condition is applied At the liquid-gas interface, z=h(x,y,t), appropriate boundary condi- 
tions are applied. They are, the kinematic equation for free-surface motion, 

ht + Q t ,t = 0 , (25) 

balance of normal stress, 

= 2 HS , (26) 

balance of tangential stresses, for a=l,2, 

M 

crtjnjt* = —9 <t t t a , (27) 
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(28) 


and the heat balance, 


0,n, = —BiO 


Here, 

rh(x,y,t) 

Q: = (Qx,Qy,0) = / U ,dz 

Jo 

is the volume flow rate vector, n is the unit outward normal vector 


n t — {—h x , —h y , 1)(1 + h x + h y ) 


2n"2 


(29) 


(30) 


and t, 1 and t 2 are orthogonal unit tangent vectors to the interface: 

<. i = (o,i,t,)(i+vr i . (3i) 

t , 2 = n,tk l t,k, , (32) 

H is the mean curvature of the free-surface 

2 H = [h xx (l + h y 2 ) — 2 h x h y h xy + h yy (l + h x 2 )](l + h x 2 + h y 2 ) 5 , (33) 

Mis the Marangom number, Bi=(hd 0 /k ) is the Biot number, S=(aodo/[3pu 2 ]) is the surface 
tension number. 

The surface-tension enters the dynamics of the problem through the force balance at 
the free-surface. The normal stress jump at the interface is balanced by surface tension 
times twice the mean curvature of the interface. In the absence of viscosity this is called 
Laplace equation which states that the pressure is larger on the concave side of the inter- 
face by the amount 2 HS. Thermocapillary is introduced through surface tangential stress 
balance by the dependence of a on T. This can either alter the capillary pressure jump at a 
particular location or introduce surface flow where fluid flows from hot end to cold end (for 
7 >0). Since the bulk of the fluid is viscous, this will also be dragged along. This is known 
as thermocapillary effect and the instabilities of this type can drive its own instability and 
does not need any external influence. Biot number Bi determines the amount of heat loss 
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at the interface Thermocapillary flow does not occur if Bi = 0 or oo. Bi = 0 corresponds 
to insulated free-surface and the interface obtains the plate temperature. On the other 
hand Bi=c o corresponds to highly conductive fluid and the interface obtains temperature 
of ambient gas. 

The analysis is done on a three-dimensional box whose base dimensions are 2% /k\ and 
27 t/& 2 where k\ and k 2 are the wavenumbers of the imposed disturbance in the streamwise 
and spanwise direction respectively. Consequently, periodic boundary condition is imposed 
for u, p and 6 m these directions. In this problem we need to satisfy three hydrodynamic 
boundary conditions at the free-surface viz , normal stress balance, tangential stress balance 
and the kinematic equation. We incorporate the normal stress balance directly into the 
momentum equations and apply the tangential stress balance as a natural boundary. Once 
the new field variables u, p and 0 are calculated, the free-surface height is updated by 
solving the kinematic equation. This procedure is called kinematic iteration. 

3. NUMERICAL SCHEME 

The governing equations (21) to (23) along with the boundary conditions are solved 
using Semi Implicit /Exp licit Finite Element Method. In this method, the viscous and 
pressure terms of the governing equations are treated implicitly and the nonlinear convective 
terms and the kinematic equation are solved explicitly. The splitting admits the use of fast 
iterative solvers and helps to minimize storage requirements. This fractional step scheme, 
based on Helmholtz decomposition theorem proposed initially by Chorin (1968) in finite- 
difference context, is best suited for the time dependent problem like the present study. 
Starting with an initial free-surface profile, the first step is to compute an intermediate 
velocity field (u” +1 ) by omitting the pressure term from the momentum equation. Viscous 
terms are treated implicitly and the convective terms are treated explicitly. Second-order 
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Adams-Bashforth Scheme and implicit Euler method are used 


' u? +1 -u, n ' 
At 


+ ~ + ~( u j “ Sj)" V/ 1 + F « n (34) 


Here, i=j=l,2,3 and the contributions from the gravity terms are included in the load vector 
F t . Since the diffusion terms are treated implicitly, these terms do not pose any restriction 
on the stability of the scheme. The time step is chosen such that the CFL condition is 
always satisfied, i.e. 


A t<C 


Ax 


mm 


Ay A z 


u -f V Gh + 


u; 


(35) 


where C is the Courant number and C < 1 for a stable scheme. The terms inside the 
square root represent the contributions from gravity and capillary waves respectively. 

The next step consists of calculating pressure from the intermediate velocity field. 
This is accomplished by projecting u" +1 into a divergence-free space. The resulting pressure 
Poisson equation is solved by satisfying the normal stress balance at the interface. 


u n+1 

„ n+l _ 

P ' u ~ A t 


(36) 


After pressure calculation, the final velocity field u, n+1 is computed by adding suitable 
contribution of the pressure field to u” +1 


u 


n+l 


— U 


7Z-f- 1 


At 


= -P, 


n+l 


(37) 


After the final velocity field is computed, the temperature field is calculated by solving 
the energy equation in similar fashion. In this one step calculation, the convective terms 
are treated explicitly using second-order Adams-Bashforth scheme and the diffusion terms 
are treated implicitly using Euler backward method 

( r+ ^ 7 r ) + *[ ;r = -fK - + 1(^ - (38) 
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The next step is to locate the free-surface h(x,y,t) that is not known a priori. Therefore, 
at every time-step, we solve the kinematic equation (Eq.25) by assuming the interface as 
a material surface i e Q^=0. The limitation of this procedure is that h(x,y,t) needs to be 
a single-valued function of x and y We solve Eq.25 using Fourier Spectral method. The 
advantage of using this method over fimte-difference/finite element schemes is that this 
method can preserve the symmetry of the geometry over a very long period of evolution. 
A specific example is heated thin film on a horizontal solid substrate where the nonlinear 
flow develops isotropically. The position of the interface is calculated explicitly as: 


'h n+1 - h n ' 
A t , 


+ Q,,, n = 0 


(39) 


Using h n+1 , new location x, n+1 of the grid points and grid velocities g, n+1 are calculated. 
It is assumed that mesh points are resting on vertical spines and are allowed to move only 
up or down depending on the local interface height. Thus the ability of moving the nodes 
in the manner we desire is due to solving the governing equations in an ALE frame of 
reference. Since the grid points are allowed to move only parallel to z-axis 


g; 


n+1 


' z n+l _ z n 


At 


(40) 


and g 2; Tl+1 and g y n+1 are always zero. The above procedure is repeated till the film ruptures 
or the free-surface equilibrates to the desired wave form 


4. RESULTS AND DISCUSSIONS 


4-1 Rupture Dynamics 

We first integrate two-dimensional governing equations to study spontaneous rupture. 
A simple-harmonic disturbance of the form 

/i(x,0) = 1 + 01 cos(kx) (41) 
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is imposed on the free- surface and its evolution is studied in time (temporal stability analy- 
sis). Here, k is the wavenumber of the imposed perturbation. The analysis will be done on 
one wavelength periodic domain. In this way, we allow the disturbances of wavelength (A, 
A/2, A/4, ...) to grow and interact nonlinearly. This is called as study of super-harmonic 
instability. Spatial stability analysis, examining growth/decay of a disturbance in space, 
and temporal stability analysis can be converted into one another through the Gaster 
transformation. We set k=k m ,k m / 2 and k m /\ where k m is the wavenumber for maximum 
growth-rate and is defined as k c /\/ 2. k c is the cut-off wavenumber and is expressed as (Joo 
et al. 1991): 



Computation is stopped at any moment when the local film thickness becomes less than 1% 
of initial mean thickness (/i(:r,t/,f)<0.01) and rupture is assumed at that spot. Beyond this 
point intermolecular forces, neglected in our formulation, become significant. However, in 
the absence of ionic molecules m the liquid, only van der Waals force of molecular attraction 
will be significant. This force is destabilizing and the film will ultimately rupture soon. 

The computational domain is discretized into non-overlapping three node linear trian- 
gular elements. Grid convergence study is conducted for different cases simulated so that 
we resolve even small scale structures of the flow. In Fig.3 results from one such study 
is shown for (7=1, 5=100, P=7.02, M— 35.1, Bi= 1.0, /?=0 and k=k m . In these figures, 
snapshot of the interface is shown at different time levels until the film breaks. In (a) we 
have used 32 modes (33 grid points in the x direction and 11 grid points in the y direction) 
to solve the kinematic equation. Similarly, we have used 64 and 128 modes respectively in 
(b) and (c). These figures clearly prove that a grid size of 65 mesh points in the x direction 
and 11 grid points in the y direction is sufficient. Subsequently, twice or four times this 
mesh size is used respectively when the wavenumber is k m /2 or k m / 4. As it is seen, the 
spectral method has the advantage that even for coarse grid (33x11), it is moderately accu- 
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(a) (b) (c) 



0 50 0 50 0 50 

XXX 

Figure 3: Grid independent study is shown for simulation are (7=1, 5=100, 

P=1 . 02, M=35.1,5i=l, £=0.0169 and f3=0. The mesh sizes are (a) 33 x 11 (b) 65 x 11 
(c) 129 x 11 


rate enough to preserve the symmetry of the flow geometry. Every problem is highly grid 
dependent and for different cases studied, different mesh sizes are used to ensure proper 
resolution of the flow field. 

When the film is horizontal, surface-wave instability does not present However, sur- 
face tension varies along the free-surface due to the presence of temperature gradient and 
sets up the Marangoni convection. If the wavelength of the imposed disturbance is suffi- 
ciently small and the film is moderately thick, the surface-tension and hydrostatic pressure 
will stabilize the flow. However, in very thm films, the interface continues to evolve for 
long surface waves and ultimately ruptures. One such case is shown in Fig 4. The pa- 
rameters are G=1,5=100,P=7.02,M=106.2, and Bi= 0.1 The snapshot of the interface is 
shown in (a) to (c) for k=k m ,k m /2 and k m /4 respectively. The shape of the free-surface is 
approximated by a finite Fourier series 

N 

h(x,t) = a n e tknx + c.c (43) 

n=—N 

and the growth of the first four harmonic modes is shown m (d) to (f) for these cases. 
Initially, the harmonics grow exponentially as predicted by the linear theory The energy 
is confined to the fundamental mode and the free-surface maintains its simple-harmonic 
configuration. As time progresses, thermocapillary becomes significant and thinning of the 
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Figure 4: <3=1, 5=100, P=7.02,M=106.18,£h=0.1,& m =0.0677 and f3= 0. In figures (a),(b) 
and (c), evolution of the free- surface at the intervals of 100,100, and 200 viscous time units 
is shown when the disturbance wavenumber is k m , k m / 2 and k m / 4 respectively and the 
corresponding growth in the harmonic modes is shown m (d),(e) and (f). Free surface 
profile at the point of rupture is also shown along with the rupture time. These results are 
obtained from finite element simulation. 


19 





film continues at the trough. The amplitude of the wave grows and the thin layer of the 
fluid flattens due to the proximity of the solid wall. The two edges of this flat layer has large 
slope and very high curvature and are rapidly drawn downwards by the capillary pressure. 
The pressure of the liquid is lower near the edges than the flat region. This causes the liquid 
to drain outwards. Meanwhile thermocapillary induces large velocities toward the plate. 
Consequently, the liquid trapped inside the flat region moves upwards to conserve mass 
This results in a characteristic bulge at the center. Now the energy is no longer confined 
to the fundamental mode but has spread to its harmonics as seen in Fig. 4d to 4f. The 
edges continue to bulge downwards and grow isotropically The growing fingers ultimately 
touch the plate (/i(x,t)<0.01) and breaks the film. When the length of the domain is longer 
(k=k m f 2 and k m / 4), we can notice several characteristic fingers Finite element mesh at 
the point of rupture is shown in Fig.5 for these cases 

In Fig. 6, the evolution of the interface for the same case from the numerical compu- 
tation of the long-wave evolution equation (Joo et al. 1991) is presented. As the harmonics 
grow by nonlinear exchange of energy, the Fourier spectrum (Fig. 6d to 6f) broadens to 
include modes outside long- wave theory. Also, the slope of the free-surface increases rapidly 
and this violates basic assumptions used in long-wave theory formulation Therefore, this 
method fails to follow the dynamics up to the point of rupture However, both full-scale and 
long-wave computations predict formation of fingers before rupture and confirm that the 
“fingering” process is not an artifact of long- wave theory but indeed an actual phenomenon. 

Thermocapillary instability is absent if Bi = 0 or Bi—oo. When Bi= 0, the interface 
temperature (T,) is same as the temperature of the bottom plate (T w ) and when Bi= oo, T, is 
same as the ambient temperature (T s ). Further, Eq.42 shows that Bi= 1 is the critical value. 
Consequently, when we increase the Biot number to 1, we increase the thermocapillary and 
accelerate rupturing process. Knshnamoorthy et al. (1995) have already confirmed this by 
the full-scale direct numerical simulation 
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Figure 5: Typical finite element mesh at the point of rupture. The parameters used in 
the simulation are G'=l,5=100,P=7.02,M=106.2,Bz=0.1, fc m =0.0677 and (3=0. (a) k m (b) 

km/ 2 (c) km/ A. 
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(a) (d) 



Figure 6 G=\, 5=100, P=7. 02, M=106. 18, Bi=0.1,k m =0 . 0677 and 0—0. In figures (a),(b) 

and (c), evolution of the free-surface at the intervals of 100,100, and 200 viscous time units 
is shown when the disturbance wavenumber is k m , k m / 2 and k m / 4 respectively and the 
corresponding growth in the harmonic modes is shown in (d),(e) and (f). These results are 
obtained from the spectral calculation of long-wave evolution equation (Joo et al. 1991) 
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When the plate is tilted, i.e., f3 > 0, the effect of the mean flow is present. Joo et 
al. (1991) have showed that when there is no surface-wave instability and thermocapillary 
is weak, the flow equilibrates after initial instability. When the surface-wave instability is 
present, thermocapillary enhanced this instability and promotes wave breaking. When the 
thermocapillary is strong, the disturbance grows continuously and results in rupture. In 
Fig. 7, thin film flowing down on a vertical plate is shown when (7=5, 5=100, Bi=Q.l,P=7 
and M= 70. In this case there is no hydrostatic stabilization and the flow is driven by gravity. 
Initially, the harmonics grow exponentially as per linear theory. As time progresses, there 
is a nonlinear exchange of energy among the modes and since thermocapillary is weak, the 
surface- wave instability causes the film to saturate to a permanent wave form. However, for 
ft=fc m /4, the film never saturates and there is a continuous exchange of energy among the 
modes. This also indicates the influence of wavenumber on the nonlinear flow development. 
Evolution of the spectral coefficients shows that at t=1138.31, the mode n=±2 dominates 
and the wave has two peaks In Fig. 8, we increase the thermocapillary (BiM/P= 10) and 
reduce the surface-wave instability ((7=1). Due to the strong thermocapillary instability, 
the film never saturates, but continue to develop and breaks. 

4-2 Rivulet formation 

We solve an initial value problem by imposing a simple-harmonic disturbance of the 

form 

h(x,y,0) = 1 + 0.1 cos(kix) + 0.1 cos [k 2 y) (44) 

where k\ and k 2 are respectively the streamwise and the spanwise wavenumber of the 
disturbance such that k=yj Aq 2 + k 2 2 . We integrate three-dimensional governing equations 
and examine the evolution of the perturbation m time. The analysis will be done on one 
spatial period box whose base dimension is [2-k / kifl-K /k 2 ). The growth of subharmonics, 
that can occur in laboratory situation, is not allowed. We select k such that k < k c where 
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(a) 


(d) 




Figure 7: <5=5, 5=171, P=7,M=70,5z=0.1,& m =0. 1103 and 0=9 0. In figures (a),(b) and 

(c) , the free-surface shape is shown for the time indicated when k=k m ,k m /2 and k m / 4 
respectively. Corresponding growth in the harmonic modes during this period is shown in 

(d) ,(e) and (f). These results are obtained from finite element simulation 
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(a) 


(d) 




Figure 8: (7=1, 5=100, P=7,M=700,Bi=0.l,k m =0 . 2049 and /7=90. In figure (a), the free- 
surface shape at the point of rupture and in (b) evolution spectral coefficients are shown 
for k=k m . These results are obtained from finite element simulation. 


k c is the cut-off wavenumber and is expressed as (Joo et al. 1993): 


k c = <- 


BiM ( 1 


P +Bi 


2 o/''i2 r< 

6KJT . 9 r, 9/1 Or r, 

+ — — sin pcos 0 — —cos p 

15 3 


(45) 


We stop the computation at any moment when the local film thickness becomes less than 
0.01 and rupture is assumed at that spot. Beyond this point intermolecular forces, ne- 
glected in our formulation, become significant. As explained in § 4.1, these forces are 
destabilizing and ultimately break the film. The computational domain is discretized into 
non-overlapping four node linear tetrahedral elements such that six of these elements make 
one cubic box. 

When the film is vertical, the flow is driven by gravity and hydrostatic pressure does 
not present. Thermocapillary and surface-wave instabilities can enhance each other and 
the film flow may saturate or rupture depending on which mode of instability is dominant. 
In our simulations, the plane is inclined only in the streamwise direction. Consequently, 
there will not be any mean flow in the spanwise direction. In Fig.9a to 9e evolution of a 
vertical falling film is shown for t=150, 300, 600, 970 and 1027 respectively when (7=1, 
5=100, Bi= 1, M— 35 and P=l. In each figure, the free-surface shape, contour plot of 
interface height predicted by finite element calculation and spectral computation of long- 
wave theory are shown in order. The wavenumbers of the initial disturbance are ki=&2=0.5. 
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Initially, as the liquid drains downward, surface-wave instability is dominant and the flow 
evolves downstream and is shown in Fig. 9a. Here, the local phase speed of the layer is 
proportional to its thickness as per linear theory. As time progresses, the local thinning 
of the liquid layer persists (Fig. 9b at f=300). Now thermocapillary begins to dictate the 
growth of the liquid layer and the transverse wave is affected by the three-dimensional 
instability. In the absence of mean flow in the spanwise direction, the liquid is displaced 
laterally (Fig. 9c) by thermocapillary instability. This process is similar to the evolution of a 
heated thin film on a horizontal substrate. The “thin layer” effect causes the fingers to grow 
and a three-dimensional longitudinal pattern (rivulet) develops along the centerline of the 
stream. At f=970 (Fig.9d), all the superharmomcs are excited by this nonlinear exchange 
of energy and the long-wave theory is no longer valid beyond this point. However, using 
full-scale computation, we can integrate the governing equations all the way to rupture. 
The final state is shown m Fig.9e. This simulation confirms Joo et al.( 1995) observation 
that the longitudinal rivulets aligned with the mean flow can form only when both the 
thermocapillary and surface-wave instabilities are properly balanced and neither of these 
two instabilities alone has the tendency to develop such pattern. This simulation also 
explains a mechanism for rivulet formation from purely a stability point of view. 

In Fig. 10a to 10c, the evolution of the thin film is shown for spanwise wavenumber 
&2=0. 25 for various time levels indicated. All other parameters are identical to the previous 
case. In this case local thinning rates are smaller and so the rupture time is increased. 
However, “fingering” occurs in an early stage of the evolution and the rivulets are much 
larger. 
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Figure 9 (e) t=1027 


Figure 9: (7=1, 5=100, £h=l,M=35,.P=7,&i =0.5, & 2 =0.5, and /?=90. The free-surface shape, 
contour plots of free-surface height from finite element simulation and spectral computation 
of long-wave evolution equation (Joo et al. 1995) are shown in order for (a) t=150, (b) 
t=300, (c) t=600, (d) t=970 and (e) t=1027. At (=1027, long-wave theory is no longer 
valid and only results from full-scale computation are shown. 
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5. CONCLUDING REMARKS 


In this study, spontaneous rupture and rivulet formation of a thin film under combined 
thermocapillary and surface-wave instability is studied by solving the complete system 
of governing equations with fully nonlinear boundary conditions. Finite element method 
based on a projection scheme is used and the governing equations are solved in an Arbitrary 
Lagrangian Eulerian frame of reference. It is shown that spontaneous rupture always occurs 
by a “fingering” mechanism as predicted by long-wave theory. The growth of the fingers is 
isotropic when the film is horizontal. When we tilt the plate, surface-wave instability sets 
in. If thermocapillary is dominant, the film ruptures. On the other hand, if surface-wave 
instability is significant, the flow saturates into a steady wave form. Besides, the amount 
of heating, thickness of the film, angle of inclination and the amount of heat loss at the 
interface, the development of secondary flow shows strong dependency on the wavenumber 
of the imposed perturbation. In three-dimensional flow, when both thermocapillary and 
surface- wave instabilities are properly balanced, longitudinal rivulets aligned with the mean 
flow forms. 
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